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Preface 


The  purpose  of  this  study  is  to  compare  the  basis  function  and 
iteration  methods  for  solving  the  Fredholm  integral  equations  of  the  first 
kind.  The  equation  arises  in  the  analysis  of  data  from  underground 
nuclear  effects  simulations.  The  Fredholm  equation  is  ill-posed  and  so  its 
answer  is  dependent  upon  the  method  of  solution. 

Often  the  Fredhom  equation  is  solved  with  the  iteration  method  tested 
here,  but  the  basis  function  method,  tested  by  Capt  Russell  Daniel,  is 
another  productive  approach. 

There  are  two  people  that  I  would  like  to  thank  for  their  help  during 
this  study.  First  I  am  Indebted  to  my  faculty  advisor,  LCDR  Kirk  A. 
Mathews,  for  his  help  and  patience  during  this  study.  The  second  person 
Is  Capt  Russell  Daniel.  His  program  and  thesis  were  essential  to  my 


understanding  of  the  basis  function  method. 
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Abstract 


>  The  purpose  of  this  study  is  to  compare  methods  for  solving  the  Fred¬ 
holm  integral  equation  of  the  first  kind.  The  Fredholm  equation  has  sev¬ 
eral  piactlcal  applications  including  geology,  superconductivity,  and 
aerodynamics.  Of  specific  interest  is  its  application  to  determining 
radiation  spectra  using  data  from  underground  nuclear  effects  simulations. 

The  two  basic  solution  methods  studied  were  the  basis  function  and 
the  iteration  methods.  The  basis  function  method  is  a  representation  of 
the  unfolded  spectrum  by  a  series  of  Planckian  or  cubic  spline  functions. 
The  iteration  method  scales  the  unfolded  spectrum  so  that  W «  weighted 
integral  over  a  given  interval  matches  that  of  the  actual  spectrum. 

Both  basis  function  methods  produced  excellent  results  when  the 
actual  spectrum  was  a  sum  of  it’s  basis  functions.  The  cubic  spline  method 
produced  unfolded  spectra  which  were  good  approximations  for  discontinuous 
actual  spectra.  However,  there  was  a  significant  dropoff  of  the  spectrum 
for  the  cubic  spline  for  higher  energies. 

The  iteration  method  produced  accurate  approximations  for  actual 
spectra  that  were  both  basis  function  and  discontinuous  spectra.  There 
were  two  problems  with  this  method:  the  unfolded  spectra  were  discontinu¬ 
ous  at  the  discontinuities  of  the  weighting  function  and  noisy  data  some¬ 
times  produced  large  discontinuities  in  the  unfolded  spectra.  v 

For  the  four  test  cases  studied,  the  iteration  method  proved  to  be  the 

best  because  of  it's  ability  to  unfold  a  variety  of  spectra  (in  contrast  to 
the  Planckian  basis  function  method)  and  it's  good  behavior  at  higher 
energies  for  noisy  data  (in  contrast  to  the  cubic  spline  method). 


v 


1.  Introduction 


Background 

The  Fredholm  Integral  equation  of  the  first  kind  has  widespread 
application  in  areas  such  as  acoustics,  optics  (3:1-1),  superconductivity, 
geophysics,  and  aerodynamics  (4:261-263).  A  particular  application  of  the 
Fredholm  equation  is  to  determine  the  radiation  spectra  of  simulated 
nuclear  weapons  effects.  The  terminology  in  this  study  originates  from  this 
application. 

The  Fredholm  Integral  equation  of  the  first  kind  has  the  following 

form: 


/?(€,Tl)s(Tl)dTl 


(1) 


where 


K(e)  =  detector  output  signal 
£(€.n)  =  detector  response  function 
s(n)dq  =  radiation  spectrum 

In  this  application,  the  practical  aspects  of  data  collection  force 
several  limitations  on  the  detectors.  The  limitations  occur  because  the 
number  of  detectors  Is  finite,  the  detectors  have  a  finite  resolution,  E/£E, 
and  recording,  transmitting,  and  calibrating  errors  affect  the  data. 
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Recording  and  transmitting  errors  are  reduced  by  using  two  types  of 
detectors,  "open"  (fluorescer  system)  and  "closed"  (fllter/fluoreacer  system) 
response  detectors.  Theoretical  modeling  and  calibration  measurements 
determine  the  response  functions  of  the  detectors. 

Problem 

The  Fredholm  equation  Is  ill- posed;  therefore,  there  are  an  Infinite 
number  of  spectra,  s(q)drv  which  would  have  the  same  detector  output, 
K(e)  The  three  methods  of  unfolding  used  in  this  study  each  produce  an 
approximate  solution  to  the  Fredholm  equation.  The  main  purpose  of  this 
study  Is  to  compare  how  well  the  approximate  spectra,  from  these  methods, 
reproduce  the  original  spectrum.  To  do  this  comparison,  all  three  methods 
of  solution  were  Incorporated  Into  a  computer  program. 

Scope 

Three  approximate  methods  for  solving  Equation  (1)  will  be  used  dur¬ 
ing  this  study:  the  basis  function  method,  the  Iteration  method,  and  the 
modified  Iteration  method. 

The  measured  signals  were  calculated  assuming  the  actual  spectrum, 
which  In  practice  would  not  be  known.  The  unfolding  methods  produce 
unfolded  spectra.  These  unfolded  spectra  are  then  compared  to  the  actual 
spectrum. 

There  are  four  actual  spectra  that  will  be  used  in  this  study  to  com¬ 
pare  these  methods:  a  twc- temperature  Plancklan  spectrum  with  a  spike 
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added  for  reasonable  detail,  a  Heaviside  function,  a  spike  on  a  flat 
spectrum,  and  a  triangular  spike  on  a  flat  spectrum. 

The  first  spectrum  represents  a  simplified  version  of  expected  nuclear 
effects  radiation;  the  basic  spectrum  being  composed  of  multiple  Plancklans 
with  the  spike  added  as  a  reasonable  detail.  The  latter  three  spectra  are 
tested  to  see  how  the  unfolding  methods  deal  with  different  discontinuous 
behavior. 

The  tests  were  conducted  first  on  unfolded  spectra  with  no  noise. 

This  indicates  how  well  the  methods  intrinsically  work.  Then  the  methods 
were  applied  to  the  same  spectra  in  a  more  realistic  environment  of  exper¬ 
imental  uncertainties. 

Assumptions 

In  order  to  compare  with  previous  work  by  Daniel  (3)  similar  assump¬ 
tions  are  made.  They  are: 

1.  There  are  20  detectors.  13  of  these  are  "closed"  and  7  are 
"open"  response  detectors. 

2.  The  resolving  power  of  the  detectors,  E/A  E,  is  approximately 

1.6. 

3.  The  response  functions  are  exact.  Their  exact  form  is  given 
in  Section  II. 
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Also,  several  tacit  assumptions  are  made  to  restrict  the  Infinite  number  of 
solutions  to  the  unique  unfolded  spectrum  given  by  the  methods  studied. 
These  depend  upon  the  method  used  and  are  discussed  In  Section  II. 

General  Approach 

The  actual  spectrum,  which  Is  In  practice  unknown,  is  assumed  In 
order  to  calculate  the  measured  signal  by  an  approximate  solution  of  Equa¬ 
tion  (1).  One  of  the  three  unfolding  methods  is  then  used  to  determine  the 
unfolded  spectrum,  given  the  measured  signal  and  the  detector  response 
function. 

For  the  basis-function  method,  a  x*  te8t  statistic  Is  used  to  stop  the 
method.  The  Iteration  methods  are,  In  practice,  stopped  when  the  spectrum 
Is  Judged  reasonable.  If  the  method  Is  stopped  too  soon,  the  effects  of  the 
smoother  have  not  yet  had  a  chance  to  work.  If  the  method  Is  stopped  too 
late,  all  the  details  of  the  spectrum  have  been  washed  out. 

The  unfolded  spectra  will  then  be  compared  to  the  original  spectrum, 
the  actual  spectrum,  and  the  spectra  of  the  other  methods. 

Sequence  of  Presentation 

The  presentation  of  the  rest  of  the  study  is  as  follows.  Chapter  II, 

"Theoretical  Discussion",  details  the  theory  of  the  unfolding  techniques 
used.  Chapter  III,  "Computer  Program  Implementation",  describes  the  pro¬ 
gram  which  Implements  the  unfolding  techniques.  Chapter  IV,  "Results  and 
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Discussion",  reviews  the  results  of  the  unfolding  techniques.  Chapter  V, 
"Conclusion  and  Summary",  summarizes  the  results  and  presents  the  conclu¬ 
sions. 
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II.  Theoretical  Development 


Introduction 

In  this  chapter  the  theoretical  background  of  several  unfolding  tech¬ 
niques  are  discussed.  The  following  topics  are  addressed:  the  detector 
response  function,  the  definition  of  the  spectra,  the  definition  of  the 
signals,  the  definition  of  unfolding  methods,  the  flux  non-negativity  con¬ 
straint,  the  formulation  of  the  x2  test  statistic,  the  minimization  of  the  x2 
statistic,  and  the  normal  distribution  for  simulating  measurement  error. 

Linear  Fredholm  Integral  Equation  of  the  First  Kind 

The  detectors,  used  to  measure  spectra  of  underground  nuclear  simu¬ 
lations,  are  designed  to  have  an  output  proportional  to  the  spectrum  in  the 
energy  range  E  to  E+&E 


yl(£)A£=/?l(£-)s(£)A£  (2) 


where 


y,(t]af  =  is  the  Incremental  output  of  detector  i  due  to  the 
spectrum  in  the  energy  range  E  to 

R,(E )  =  response  of  detector  l  to  the  part  of  the  spectrum  at 
the  energy  E 

s(F)Ai  =  the  amount  of  the  spectrum  that  is  in  the  energy 
range  E  to  f  t  £  £ 

The  output  signal  of  detector  i  due  to  the  total  spectrum  would  be: 
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yr  f  y XE)d£-  f  Rt{E)s{E)dE 

J  0  Jo 


(3) 


where 

y.=  the  output  signal  of  detector  i  in  response  to  the  entire 
energy  spectrum. 

This  is  a  linear  Fredholm  integral  equation  of  the  first  kind. 

Detector  Response  Function 

The  detector  response  function,  /?,(£),  is  the  normalized  output  of 

detector,  1,  due  to  the  portion  of  the  spectrum  at  the  energy  E.  These 
response  functions  are  calibrated  to  a  known  spectrum  and  are  assumed  for 
this  study  to  be  the  same  for  the  unknown  source  being  measured. 

In  a  realistic  experiment,  multiple  detectors  (response  functions)  are 
used  simultaneously  and  their  results  are  compared  in  an  effort  to  reduce 
measurement  errors.  There  are  two  types  of  response  functions  used  in 
this  study:  "open"  and  "closed".  The  "open"  detector  system  is  a  fluorescer 
detector  that  has  a  zero  response  for  energies  less  than  the  k-edge  of  the 
fluorescer  and  an  exponentially-decreasing  response  for  higher  energies  as 
shown  in  Figure  I. 

The  equation  for  the  "open"  response  function  (3:3-3)  is  given  by: 


R°pt\E) 


0  for  £<£? 


for  £££t° 


(4) 


where 


=  tjxe  response  function  for  the  "open"  detector  i 

£?  =  the  k-edge  of  the  fluorescer  for  detector  i 
E  =  energy 

Figure  I  is  an  example  of  an  "open"  response  function.  The  sensitivity  is 
the  amplitude  of  the  response  function  of  "open"  detector  15.  The  energy 
axis  is  in  units  of  the  k-edge  of  the  fluorescer  of  detector  1,  EO.  All 
energy  units  in  this  study  are  given  in  units  of  EO. 


8 


The  "closed"  detector  system  is  a  filter-fluorscer  detector  which 
shapes  the  "open"  detector-type  response  into  a  more  pulse-like  envelope. 
This  envelope,  or  Inband  response,  is  between  the  k-edges  of  the  filter  and 
the  fluorescer.  Both  detector  systems  studied  are  asymmetric. 

The  equation  for  the  closed  response  function  (3:3-3)  is  given  by: 
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(5) 


0  for  E<E°l 


Rc,‘tnd(E ) 


-(^){I“exp[“2(y)3  }exp["°-25(y)3]  for  E'*E<E' 
(f){'-'*p[-2(r)  ‘^(r)  ]  for  e>£! 


Figure  II  is  an  example  of  a  closed  response  function.  The  sensitivity  is 
the  amplitude  of  the  response  function  of  "open"  detector  15.  The  energy 
axis  is  in  units  of  the  k-edge  of  the  fluorescer  of  detector  1,  EO. 
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Both  the  "open"  and  "closed"  response  functions  discussed  above  are 


Ideal.  In  the  presence  of  calibration  and  measurement  errors,  the  real 
response  function  is  perturbed  In  the  following  way: 

(6) 


R,(E)-R,(E)  +  ARi(E) 


where 


£.(£)  =  the  actual  response  function 

A/?,(F)  =  the  change  in  the  response  function  due  to  calibration 
and  measurement  errors.  This  change  could  be  negative  giving 
a  smaller  value  for  the  response  function. 

Definition  of  Spectra 

Three  spectra  are  used  in  this  study:  the  predicted  spectrum,  the 
actual  spectrum,  and  the  unfolded  spectrum.  The  predicted  spectrum, 
Sp(£),  is  the  spectrum  expected  from  theoretical  calculations.  This  spec¬ 
trum  is  known. 

The  actual  spectrum,  S„(£),  is  the  spectrum  to  which  the  detectors 

respond.  In  the  field  tests,  the  actual  spectrum  is  unknown  and  to  some 

extent,  unknowable.  For  the  purpose  of  comparing  the  performance  of  the 

unfolding  methods,  s„(f)  is  chosen.  The  ability  of  the  unfolding  methods 
to  recover  sa(f)  is  the  measure  of  the  method's  suitability. 
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The  unfolded  spectrum,  SU(G),  is  calculated  from  detector  output  and 
response  function  data  to  be  an  approximation  to  S«(E).  The  predicted 
spectrum  may  be  used  as  an  initlai  guess  for  the  unfolding  methods. 

Definition  of  Signals 

A  signal,  Yi  produced  by  a  spectrum  s(E)dE,  is  the  detector  output 
given  by  the  following  equation: 

r,«  /”i?l(f)s(f)d£  (7) 

Jo 

There  are  four  signals  associated  with  the  three  spectra  already 
defined:  the  ideal  signal,  the  measured  signal,  the  predicted  signal,  and  the 
unfolded  signal.  The  ideal  signal,  Yf*a\  is  the  signal  produced  by  the 
actual,  but  unknown,  spectrum  as  measured  by  an  ideal  instrument.  The 
ideal  signal  is  unknown  because  of  errors  in  calibration  and  uncertainties 
in  measurement  and  detection.  The  ideal  signal  is  given  by: 

r;d8ai-  f  Rl(E)Sa(£)d£  (8) 

Jo 

The  measured  signal,  y",0,ttr,<<t  is  the  signal  measured  by  the  detectors 

with  ail  the  calibration  errors  and  measurement  and  detection  uncertainties 
Included.  The  measured  signal  is  given  by: 
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Y  m.«asurtd  m  /*  £  (  (  £  )  5  a  (  £  )  d  £  (9) 

0 

The  predicted  signal,  Yfr*dut*d,  is  the  signal  expected  from  the  predicted 

spectrum  as  measured  by  the  Ideal  Instrument.  The  predicted  signal  Is 
given  by: 


y  predicted.  m  j  £,(£)£  {E)dE  (10) 

J  0 

Finally  the  unfolded  signal,  y?*'0lMi  is  the  signal  that  would  be  pro¬ 
duced  by  the  unfolded  spectrum  as  measured  by  the  ideal  instrument.  The 
unfolded  signal  is  given  by: 

y  unfoidtd  m  f  /?((£)5u(£)d£  (H) 

Jo 

The  signal  data  from  underground  nuclear  effects  simulations  are  fre¬ 
quently  specified  as  a  ratio  to  the  predicted  signal;  it  is  in  this  form  that 
the  x 2  test  statistic  will  be  formulated.  The  measured-to-predicted  ratio 
is  given  by: 


b 


Y 

Y 


m«asur*d 

l 

prQdictid 

L 


and  the  unfolded-to-predlcted  ratio  is  given  by: 


(12) 
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Ct- 


Y  unfoldtd 


rt 


prtdlcttd. 


(13) 


Errors  in  the  Measurement  Process 

There  are  two  basic  errors  that  occur  during  the  measurement  process: 
measurement  error  and  calibration  error.  The  measurement  error  is  the 
uncertainty  in  the  measurement  of  the  detector  output.  This  uncertainty 
includes  errors  introduced  in  the  transmission  of  the  signal,  the  recording 
of  the  signal,  and  the  reading  of  the  signal. 

The  calibration  error  is  the  uncertainty  in  the  determination  of  the 
detector  response  function.  Experiments  are  conducted  to  determine  the 
detector  response  function  by  using  known  x-ray  sources.  These  calibra¬ 
tions  are  subject  to  errors  in  the  measurement  process. 

The  experimenter  accounts  for  these  errors  by  the  specifying  a 
measurement  standard  deviation,  a,  for  each  of  the  detectors.  These  devi¬ 
ations  are  based  on  past  experiments  and  detector  calibrations. 

Normal  Distribution  for  Simulating  Measurement  Error 

The  measurement  errors  are  assumed,  in  this  study,  to  be  normally 
distributed.  They  are  simulated  by  random  number  sampling  of  a  Gaussian 
probability  distribution  function,  which  is  approximated  by  a  sum  of  con¬ 
tinuous  uniform  distributions  l.e.,  a  sum  of  random  numbers.  Let  X  be  a 
random  variable  which  is  defined  by: 
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(14) 
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X  -  I  xn-6 

n  -  1 


where 


X„€i/[0,  l] 

U[0 ,  l  ]  is  a  uniform  distribution  on  the  interval  |0,l]. 


The  mean  for  X  is  (8:565): 

[ix~  E[X]»  E[12U[0 ,  1  ]  —  6 ]  =  \  2E[U[ 0,  l]]-6-0  (15) 


Since  o2  for  a  uniform  distribution  of  the  interval  [0, 1 }  is  given  by  the 
following  (8:570): 


0 


1 

12 


(16) 


and  the  standard  deviation  is  (8:566): 


0  Y  *  0 


i.*-  120, 


X  v  124/(0.  I  J-6  1  ^  v  4/(0.  1  ] 


* 12'  Ti"  1 


(17) 


then 


X*C(ja-0;o-  1) 


(18) 


which  is  the  standard  normal  distribution. 
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Definition  of  Unfolding  Methods 

Given  that  the  measured  signal  (y »*«**«<!  )(  predicted  spectrum  (SP(F)  ), 

and  the  response  functions  for  the  detectors  (£,(£)  )  are  known,  then  the 
problem  is  to  solve  for  the  actual  spectrum,  Sa(F),  in  the  system  of  equa¬ 
tions: 


r  Ri(E)Sa(E)dE  (19) 

J  0 

However  this  problem  is  ill-posed  (6:175),  in  that  there  are  an  infinite 
number  of  spectra,  S0(F),  which  will  satisfy  this  equation.  So  additional 
assumptions  are  imposed  to  determine  a  unique  solution,  S„(F),  such  that 
su(F)-sa(F).  The  unfolded  signals  should  equal,  or  at  least  approximate, 
the  measured  signals.  The  closeness  of  the  approximation  depends  on  the 
estimation  of  the  accuracy  of  the  measurement  and  is  discussed  below.  So, 

yrnta^na-  f  £  (  )  5  a  (  £  )  d  £  *  f  R  .  (  £  )  5  ^  (  £  )  d  £  (20) 

J  0  Jo 

or,  in  order  to  solve  numerically 

Y^a^  =  klRikS  AEk-kTRL,kSakAEk  (21) 

k-0  K  k-0 

The  other  constraints  are  dependent  upon  the  method  of  solving  for  SU(E). 
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Buis  Function  Method 


There  are  three  methods,  examined  in  this  study,  to  compute  a  "rea¬ 
sonable"  spectrum  for  Su(E):  basis  function(s)  expansion,  iteration  of  a  bet¬ 
ter  fit  and  smoother  operation,  and  iteration  of  a  better  fit  using  a 
smoothed  response  function.  The  basis  function  expansion,  with  Planklan 
and  cubic  spline  basis  functions,  was  studied  by  Daniel  (3:1-4).  The  idea 
is  to  approximate  $S{E)  as  a  sum  of  basis  function(s),  so  S„(E)  is  chosen 
to  satisfy: 


Sa(£) 


max 


I  cl.B  ,  k 

/-i  '  “ 


(22) 


where 


j  =  summation  index  over  basis  functions 
B  B  .(E  -,q)  is  the  basis  function  j 
q  =  the  defining  parameters  for  Bj 

The  parameters,  q,  define  the  basis  function  and  are  dependent  upon 

the  type  of  basis  function.  For  Plankian  basis  functions  these  parameters 
reduce  to  a  single  value,  the  temperature  of  the  Plankian.  However,  for 
cubic  splines,  the  parameters  include  the  locations  of  the  knots  and  the 
temperature  of  the  Plankian  fit  between  the  next-to-last  and  last  knots. 
Substituting  Equation  (22)  into  Equation  (21)  gives, 
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Y 


unfolded 

L 


l  Rltk  I  [a,B j  k]&E k 

fc» 0  /-I 


(23) 


or  rearranging, 


Y\ 


unfolded.  _  j* 

/-» 


i.  / 


(24) 


where 


6,.,-  r«,.,5(.,A£k  (25) 

k-0 

In  this  method,  a  guess  is  made  for  the  initial  unfolded  spectrum  by 
specifying  the  parameters  (a;  and  g)  and  number  of  basis  functions 
(;„„).  As  a  consequence  B,.k  can  be  determined  and  since  R,  k  and  A£*  are 
known  then  fi,,(can  be  calculated.  Finally  / can  be  evaluated  by 
using  Equation  (24).  . 

This  initial  guess  is  modified,  by  changing  a,  and/or  q  to  minimize  the 
X2  difference  between  and  yp«*‘“r*d.  The  method  is  terminated  when 

the  x2  test  statistic  of  the  difference  is  sufficiently  small.  I.e., 

2  r  y  measurement  y  unfolded  \  ^  2  (  O 

X  V'  t  •  •  L  )  ^  X  acceptable  ) 

The  formulation  of  the  v2  test  statistic  is  discussed  below. 
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Iteration  Method 


The  second  procedure,  to  compute  a  reasonable  Su{£),  is  the  iterative 

method.  The  iterative  method  is  the  repeated  application  of  two  spectrum 
changing  operations:  the  "better  fit”  and  "smoother".  The  "better  fit"  algo¬ 
rithm  forces  the  unfolded  signal,  for  the  energy  Interval  E  to  E  +  d£,  to 
equal  that  of  the  measured  spectrum: 


R ,(E  )sn"‘"‘r‘d(E)d E  ■  )df  (27) 

This  is  accomplished  by  scaling  the  guessed  spectrum,  for  that  energy 
interval,  by  the  ratio  of  the  measured  to  unfolded  signal. 


R,(E)s{r'‘](E)dE  = 


y  measured 

~^RXE)sl'\E)dE 


(28) 


where. 


2  (F)df  is  better-fit  of  interation  n  of  the  unfolded 
spectrum. 

i'u,'\EjdE  is  iteration  n  of  the  unfolded  spectrum. 

Equation  (28)  can  be  rewritten  by  using  Equations  (12)  and  (13): 

RXE)s'r'l](£)dE  =  SitRt(E)sl"(E)dE  (29) 
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where 


p.- 


f'i t(£)Sa(£)d£ 

J  0 

T  Rt(E)Su(E)dE 

Jo 


>'r 

vr 


Y?/Y> 

VTTyJ 


b, 

c( 


(30) 


For  the  underground  nuclear  effects  simulations,  the  detectors  are 
designed  to  have  overlapping  energies  to  reduce  calibration  and  measure¬ 
ment  errors.  So.  the  total  unfolded  signal  for  an  energy  interval  for  all 
the  detectors  is  given  by 


I  Rt(E)S[  (E)dE 

l-  i 


lBIX 

I  3,ff,(£)S‘"’(£)cl£ 


(31) 


but  sir?J(£)d£  and  S^(F)dF 


are  not  dependent  upon  the  detector  index 


so. 


1  max 

I  Rt(E) 

_  i- 1 

finally,  the  effect  of  the  "better  fit"  on  the  last  spectrum  is  given  by: 


\  n  *  - 


SI  (E)dE  = 


I  M<(£) 

(- 1 


S(ua\E)dE  (32) 


sl"'^(£)d£- 


l-  1 


s 


('l) 

u 


(E)dE 


I  /?,(£) 


(33) 
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The  detector  response  functions,  have  discontinuities  about  the 

k-edges  of  its  filter  and  flourescer,  FIGURES  I  and  II.  During  the  "better 
fit"  algorithm,  the  scaling  factor  in  Equation  (33)  can  be  significantly  dif¬ 
ferent  even  between  adjacent  detectors. 

FIGURE  III  gives  an  Illustration  of  such  a  case,  the  obvious 
discontinuities  are  at  16  and  24  Eo.  Unrealistic  features  that  are  a  func¬ 
tion  of  the  method,  such  as  the  discontinuous  slope,  are  called  an  artifact 
of  the  method. 
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The  second  spectrum-modifying  operation  of  the  iteration  method,  the 
smoother",  is  designed  to  spread-out  the  artifact  of  the  "better-fit"  algo¬ 


rithm  so  that  the  resulting  spectrum  is  more  realistic. 

The  averaging  of  the  artifacts  with  the  surrounding  spectrum,  is  done 
in  the  following  manner: 


S‘"*l>(F)d£- 


/; 


f(E,E')S\  *'(F')df' 


/; 


-dE 


/(F.F')dF' 


(34) 


or  written  to  be  solved  numerically, 


sr° 


dE  * 


(it-U 

I  fkkS i  *'A£Y 

* -i 


df 


5:  Af* 


(50) 


produces  a  spectrum  that  is  nearer  to  S0(F)dF  than 


{E)dE. 


There  were  several  constraints  placed  upon  the  averaging  kernel,  fk.k, 

during  this  study.  For  computational  ease,  the  kernel  was  normalized  i.e., 
for  all  E 

f  f  (E ,  E  ')dE  '  ■  1  (36) 

J  o 


or  for  discrete  energy  bins, 
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(37) 


£  fk.k  A£k  ■  1 

k-  1 


This  simplifies  Equation  (50)  to 


Q(n*n 
J  k 


dE  = 


k=« 


Z  / 


k  -i 


k.k 


l  n- 

J  * ' 


I  ^ 
z) 


A  E  i 


(38) 


The  second  constraint  Imposed  on  fkk  is  that  fkk.  is  symmetric  about 
k  l.e.,  for  constant  energy  bins: 


/  k.  k-5k  m  /  fc,  k*6k  (39) 

This  is  equivalent  to  saying  that  the  averaging  algorithm  does  not  favor 
higher  or  lower  energies.  For  energy  bins  that  are  geometrically  spaced 
l.e.,  for  any  k  the  following  is  true: 

=  (40) 


where  r  is  a  constant,  then  the  symmetry  of  fkk  is  expressed  as: 

/  k.k/rm  /  k.rk  (41) 

Finally  because  fk  k.  is  symmetric  about  k,  a  symmetric  difference  ker¬ 
nel  was  often  used  to  improve  the  computational  efficiency  of  the  "smooth¬ 
ing"  algorithm  l.e.. 
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/  k.k'  m  f  \k-k  | 


(42) 


This  is  generally  much  easier  to  compute  because  the  interval 

[~6fcn.»x of  Equation  (58)  is  much  smaller  than  the  interval  [l.fc.,x]. 

The  iteration  method  is  employed  in  the  following  way.  An  initial 
spectrum,  S9U.„(E )d£  •  Si0)(£  )d£,  is  guessed.  The  "better-fit"  algorithm, 
Equation  (49)  is  applied  to  Si0>(f)d£to  produce  s[  ^(£)d£.  The 
smoother,  Equation  (58),  is  then  applied  to  produce  S£!:(£)d£  the  first 
iteration.  At  this  point,  a  judgement  would  have  to  be  made  to  evaluate  if 
the  first  iteration  is  reasonable.  If  it  is,  then  the  first  Iteration  is  the 
answer.  If  not,  then  the  Iteration  method  is  reapplied  until  one  of  the 
spectra  is  reasonable. 

In  practice,  the  "better-fit"  and  "smoother"  algorithms  compete  against 
each  other;  the  "better-fit"  pushes  the  unfolded  spectrum  toward  the  actual 
spectrum  while  the  "smoother”  pushes  the  unfolded  spectrum  away.  This 
leads  to  the  last  assumption  for  the  iteration  method: 
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The  final  unfolded  spectrum,  which  is  a  compromise  between  the  "bet¬ 
ter  fit"  and  "smoother"  algorithms,  is  a  good  approximation  to  the  actual 
spectrum.  I.e.,  there  is  a  best  N  (which  should  appear  in  the  early 
Iterations)  such  that: 

S(uN)(t)dt  *  S a(t)dE  (44) 

The  third  method  Is  the  iteration  of  a  better  fit  using  a  response 
function  smoothed  only  around  the  k-edges.  This  method  is  related  to  the 
second  method,  with  the  "smoother"  algorithm  buried  in  the  initial  smooth¬ 
ing  of  the  response  function.  However,  in  this  case,  the  "smoother"  is 
local  and  weak  compared  to  the  "better  fit"  algorithm,  so  the  overall  effect 
is  primarily  due  to  the  "better  fit"  algorithm. 

Modified  Iteration  Method 

In  order  to  remove  artifacts  yet  not  work  against  the  "better-fit" 
algorithm,  a  third  method  is  discussed  which  is  a  modification  of  the  itera¬ 
tion  method.  The  modified  Iteration  method  works  as  follows:  the  response 
function  is  smoothed  only  around  the  k-edges  of  the  detectors  i.e., 

5  k 

Ri,k~  ^  f\k-k  \R,k  &Ek  (45) 
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and  then  normalized.  In  other  words,  the  "smoother"  has  already  been  used 
on  the  response  function  to  prevent  the  discontinuities  during  the  "better- 
fit"  algorithm.  The  "better-fit"  algorithm  now  has  the  following  modified 
form: 


"t___ 

5l’,*')(F)d£=  Fa, /?,(£) 

_  i-  I 


S‘“>(F)dF 


(46) 


Most  of  the  previous  discussion  of  the  Iteration  method  also  applies  to  the 
modified  iteration  method. 


Errors  in  the  Unfolding  Process 

There  are  two  errors  associated  with  the  unfolding  process:  approxima¬ 
tion  of  an  Integral  by  a  summmation  and  the  ill-posedness  of  the  Fredholm 
integral  equation.  Whenever  an  Integral  is  replaced  by  a  summation  in 
order  to  be  numerically  computed  there  is  an  error  associated  with  the 
replacement: 


E)s(E)d 


E 


l  max 


I  Rt 

/-  ‘ 


s/A£,  +  6(47) 


where, 


R 


i.  / 


6  Is  the  error  associated  the  approximation 


O  is  dependent  upon  the  bln  size,  A  E. 

Another  error  is  due  to  the  ill-posedness  of  the  Fredholm  Integral 
equation.  Because  the  Fredholm  Integral  equation  of  the  first  kind  is  ill- 
posed,  there  are  an  infinite  number  of  the  solutions  that  satisfy  the  inte¬ 
gral  l.e., 


f  R i(E)s  \  (E)dE  =  f  Ri(E)s2(E)dE  (48) 

Jo  Jo 

doesn't  guarantee  that  s,(e)-s2(r>  The  spectrum  determined  is  dependent 
upon  the  unfolding  method. 

The  actual  spectrum  is  not  expected  to  be  determined  by  the  any  of 
the  methods  presented  here.  The  unfolding  methods  produce  a  particular 
solution  to  the  integral  equation  which  has  a  reasonable  shape  e.g.,  a  con¬ 
tinuous  slope  everywhere.  Different  methods  can  be  used  to  produce  rea¬ 
sonable  shapes  and  some  of  the  methods  will  produce  reasonable  shapes  for 
certain  detector  signals. 

Plux  Non-negativity  Constraints 

Negative  Intensity,  for  radiation,  Is  meaningless.  However  during  the 
unfolding  process,  the  spectra  can  become  negative  at  some  energies.  So  a 
non-negativity  constraint  is  applied  such  that: 
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Formulation  of  v2  Test  Statistic 

The  test  statistic  used  in  terminating  the  basis  function  method  uses 
the  "goodness  of  fit"  test  of  the  measured-to-predicted  ratios,  bi,  to  the 
unfolded-to-predicted  ratios,  ci,  and  is  given  by: 


instruments 
2  r* 

X  ■  X 

i- 1 


(50) 


where 

o,=  the  standard  deviation  of  the  measured-to-predicted  ratio, 
bi,  based  on  estimates  of  detector  calibration  and  measurement 
uncertainties. 

When 

X2*v  (51) 

where 

v  =  number  of  degrees  of  freedom  of  the  unfolding  method 
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then  the  unfolded  and  measured  spectra  are  statistically  equivalent.  For 
the  basis  function  method,  the  number  of  degrees  of  freedom  is  the  number 
of  parameters  describing  the  basis  functlon(s)  i.e.,  a,  and  qt. 

Minimization  of  y2  Test  Statistic 


To  improve  the  fit  of  the  unfolded  spectrum  calculated  by  the  basis 
function  method,  there  is  a  systematic  method  to  minimize  x1  by  varying 
the  coefficients  and  parameters:  a,  and  q,.  The  Flecher-Powell  (F-P) 
minimization  technique  was  used  to  determine  the  approximate  direction  of 
the  minimum  for  x2  1°  the  basis  functions'  parameter  space.  The  F-P  tech¬ 
nique  takes  an  initial  point,  P,  in  the  parameter  space  (which  is  determined 
by  the  unfolded  spectrum)  and  calculates  the  direction  in  the  following 
way: 


then  the  unit  normal  direction  is  given  by, 


S  = 


S' 


I 

I 


S' 


I 

I 


(52) 


(53) 


where  H  is  Initially  a  unit  matrix  of  the  parameter  space. 

Once  the  direction,  S,  has  been  approximately  determined  by  the  F-P 
technique  then  an  appropriate  distance,  w,  in  the  parameter  space  has  to 
be  determined  to  minimize  x*-  This  is  determined  in  the  following  manner: 
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the  point,  P,  in  parameter  space  has  a  value  of  ^2(/>)and  a  slope  of 
(<J/<Ju<)[x2(P)]  along  the  direction  S.  A  smali  step,  w',  is  taken  along  S 
from  P  to  a  new  point,  P  +  ui'S.  The  value,  x2(P+ w'S),  and  the  slope 
(3/3u/)[X2(/> u/'S)]  are  determined  at  the  new  point.  With  the  two  points 
and  slopes,  a  cubic  fit  approximating  the  slice  in  parameter  space  can  be 
calculated  and  a  reasonable  distance,  w,  is  determined  for  the  closest 
approach  to  the  minimum  along  that  slice. 


For  the  next  use  of  the  F-P  technique,  the  unit  matrix  H  is  changed 
toward  the  minimum  by  the  following  formulae  (6:76): 


30 


H • H+ A- B 


(54) 


where 


rrt 

Am  T 

RtQ 

(55) 

n_HQQTH 
*  QT HQ 

(56) 

p  a  x2  (P  +  w '  s )  ~  x2  (P) 

(57) 

n  -  /  i  /  i  r  }  /  } ...  \ 2  i 

\p.w  s -yvi  UW)^  1  p 

W'  ) 

If  x2  is  still  too  large,  then  the  F-P  technique  is  reapplied  until  the  result 
Is  small  enough. 
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in.  Computer  Program  Implementation 


Introduction 

With  the  theoretical  background  developed,  we  can  now  discuss  imple¬ 
mentation  of  the  unfolding  algorithm.  In  this  section,  programming  consid¬ 
erations  are  discussed. 

Input  of  the  Test  Data 

Data  provided  by  a  nuclear  effects  test  includes:  the  predicted  spec¬ 
trum,  S „(£'),  the  calibrated  detector  response  functions.  /?,(£),  the 
measured-to-predicted  ratios,  6.  and  the  detector  standard  deviations,  a. 
These  variables  are  simulated  by  specifying  them  in  a  self-consistent  man¬ 
ner. 

#,(£)  is  assumed  to  be  exact  and  given  by  Equations  (4)  and  (5).  The 
modified  iteration  method  alters  these  response  functions  as  part  of  its 
unfolding  process,  however  the  original  response  functions  are  used  in  the 
calculation  of  the  x2  test  statistic. 

The  detectors'  standard  deviation  is  specified  as  0.15  by  compatibility 
with  Daniel's  initial  assumptions  (3:1-3).  Further  discussion  of  the  mea¬ 
surement  error  is  covered  in  section  II. 

In  reality,  the  predicted  spectrum  is  based  upon  theoretical  treatment 
of  the  physics  of  the  actual  simulation  configuration  and  results  from  past 
slmufatlons.  This  spectrum  should  be  a  very  good  approximation  of  the 
actual  spectrum,  which  is  unknown.  For  this  study,  the  predicted  specta 
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are  always  a  good  approximation  to  the  actual  spectra. 

Although  the  measured-to-predicted  ratios  could  have  been  chosen 
directly,  the  specification  of  the  actual  spectrum  and  calculating  6,  with 
Equations  (10),  (11),  and  (12)  allow  a  direct  comparison  of  Sa(£)  to  Su(£). 


Basis  Functions 


Each  of  the  spectra  can  be  specified  either  as  a  basis  function  or 
explicitly  listed  data  in  a  file.  There  are  two  types  of  basis  functions 
used  in  this  study:  Planckians  and  cubic  splines.  The  Planckian  basis 
functions,  B,.k,  as  used  in  Equation  (24)  are  written  as: 


l  "exp!  - j-  1 


(58) 


where. 


E  =  center  bin  energy  in  units  of  E 0 
T  =  temperature  of  Plankian  in  units  of  E0 
k  =  Boltzmann's  constant 

The  cubic  spline  basis  functions  are  a  composite  of  four  cubic  splines 
connecting  each  knot  to  it's  surrounding  knots.  Each  of  the  splines  are 
given  by  the  Lagranglan  function: 
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£♦.,(£)-  fl 


E  -  knot , 


i-o. i+kknotk-  knott 


(59) 


where, 


E  =  center  bin  energy,  in  units  of  E0 

knoti  =  energy  of  a  surrounding  knot,  in  units  of  £  0 

knotk  =  energy  of  the  reference  knot,  in  units  of  E0 

The  first,  second,  and  last  knots  are  located  at  0,  £0,  lOf  0  respec¬ 
tively,  with  the  weights  of  the  first  and  last  knots  always  set  to  0. 
Modifications  are  used  for  knots  near  the  edges  of  the  energy  interval. 

If  knotk  has  adjacent  knots  that  are  the  first  or  second  to  last  then, 

because  a  cubic  spline  could  not  be  calculated  to  Join  the  remaining  knot, 
they  are  joined  by  linear  curves.  The  second-to-last  and  the  last  knots 
are  joined  by  a  Planckian  for  the  same  reason. 

Iteration  Method 

The  iteration  methods  use  two  different  schemes  to  smooth  out  the 
spectrum  after  the  better-fit  algorithm,  flat-averaging  and  normal¬ 
averaging.  The  flat-averaging  method  gives  equal  weighting  for  all  energy 
bins  within  the  specifed  interval.  The  normal-averaging  method  gives  a 
normal  weighting  of  the  energy  bins  about  a  specified  point. 


Both  methods  compensate  for  intervals  too  close  to  the  end  of  the 


energy  interval  (EO  to  128EO)  by  truncating  the  long  side  of  the  averaging 
interval  so  that  it  is  symmetric  about  the  central  point  of  the  Interval. 
Otherwise  the  smoothing  function  would  bias  the  unfolded  spectrum  unreal¬ 
istically  away  from  the  ends  of  the  energy  Interval.  The  two  smoothing 
schemes  incorporate  the  possibility  of  nonconstant  spacing. 

In  order  to  reduce  the  run  time  of  the  iteration  method,  the  linear 
spacing  of  the  energy  bins  was  changed  to  a  geometric  spacing.  This 
reduced  the  time  to  calculate  the  integrals  involved  in  the  measured-to- 
predicted  and  unfolded-to-predicted  ratios.  The  overall  run  time  was 
reduced  by  a  factor  of  seven. 

The  modified  iteration  and  basis  function  methods  were  left  with  linear 
energy  bin  spacing.  This  was  done  for  the  following  reasons.  The  other 
unfolding  methods,  with  linear  spacing,  had  comparable  running  times  to  the 
geometric  energy  spacing  version  of  the  iteration  method.  The  linear  spac¬ 
ing  was  kept  as  a  check  against  the  results  of  the  iteration's  geometric 
energy  bln  spacing. 

The  stopping  point  for  the  iteration  methods  is  based  upon  the  judge¬ 
ment  of  the  person  running  the  program.  It  is  almost  an  art  to  stop  the 
program  when  it  has  nearly  converged  to  the  actual  spectrum,  but  has  not 
gone  too  far  as  to  completely  smooth  out  the  original  data.  For  the  pur¬ 
poses  of  comparison  in  this  study,  both  iteration  methods  were  stopped 

after  five  iterations.  This  is  comparable  to  the  number  of  iterations  used 
in  the  field  (7). 
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IV.  Results  and  Discussion 


Introduction 

This  Section  covers  the  cases  run  using  the  unfolding  program  dis¬ 
cussed  in  the  last  section.  The  topics  covered  are  the  following:  the  vali¬ 
dation  cases,  the  effect  of  noise  on  the  validation  cases,  the  tests  cases, 
and  the  effect  of  noise  on  the  test  cases. 

Validation  Cases 

The  validation  cases  are  to  ensure  that  the  program  is  running  pro¬ 
perly.  A  "perfect"  unfolding  technique  would  produce  an  unfolded  spectrum 
is  equal  to  the  actual  spectrum.  This  perfect  reproduction  of  the  original 
spectrum  can  not  be  required  because  of  the  111-posedness  of  the  Fredholm 
equation,  as  discussed  in  section  II.  However,  the  unfolded  spectrum 
should  be  a  close  approximation  to  the  actual  spectrum.  The  initial  guess 
for  the  unfolded  spectrum  was  the  same  for  the  iteration,  modified  Itera¬ 
tion,  and  the  Planckian  basis  function  methods. 

For  the  iteration,  modified  iteration,  and  Planckian  basis  function 
methods,  the  validation  cases  consisted  of  three  cases: 

1.  The  actual  and  guessed  spectra  are  different  Planckians. 

This  is  to  show  that  the  method  will  produce  unfolded  spectra 
that  converges  toward  the  actual  spectrum.  The  actual  spectrum 

was  T  =  3E0  and  coefficient  =  0.75.  The  guessed  spectrum 
was  T  =  2E0  and  coeffient  =  0.5. 


2.  The  actual  and  guessed  spectra  are  the  same.  This  is  to 
show  how  an  initially  correct  guess  is  changed  by  the  unfolding 
method.  The  actual  and  guessed  spectra  were  T  =  3E0  and 
coefficient  =  0.75. 

3.  The  actual  spectrum  is  a  two  temperature  Pianckian  spectrum 
and  the  guessed  spectrum  is  a  one  temperature  Plankian 
spectrum.  The  actual  spectrum  was  T1  =  3E0,  coefficient  1  = 

0.75,  T2  =  10,  and  coefficient  2  =  5.  The  guessed  spectrum 
was  T  =  10  and  coefficient  =  10. 

For  case  1,  the  modified  iteration  method  produced  an  unfolded  spec¬ 
trum  that  follows  the  actual  spectrum,  but  has  discontinuities  near  the 
k-edges  of  the  response  functions,  FIGURE  V.  This  is  characteristic  of  the 
iteration  method. 
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The  unfolded  spectrum  from  the  iteration  method  had  smoother  and  had 
less  abrupt  discontinuities  than  the  modified  iteration’s  solution  did. 
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The  Planckian  basis  function  method  gave  an  unfolded  spectrum  that 
appeared  Identical  to  the  actual  spectrum.  All  three  methods  successfully 
produced  an  unfolded  spectrum  which  a  good  approximation  to  the  actual 
spectrum,  when  the  initial  guess  was  wrong. 

For  case  2,  both  iteration  methods  and  the  Plankian  basis  function 

method  produced  an  unfolded  spectrum  that  was  nearly  identical  to  the 

actual  spectrum  (no  figure  was  produced  for  this  case  because  the  spectra 
could  not  be  distinguished).  In  other  words,  the  unfolding  methods  didn't 

distort  a  correct  guess. 
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For  case  3,  both  iteration  methods  produced  over-/under-shots  for  the 


greater  energies,  Figures  VII.  The  Planckian  method  produced  a  good 
representation  of  the  actual  spectrum,  as  expected  for  a  method  to  repre¬ 
sent  Planckian  spectra. 


FIGURE  VII:  Modified  Iteration  Method,  Case  3 


For  the  cubic  spline,  the  validation  cases  consisted  of  two  cases: 

1.  The  actual  and  guessed  cubic  splines  were  different.  The 
actual  spectrum  had  six  knots  at  the  following  locations:  0, 
EO,  5E0,  30E0,  60E0,  and  128E0.  The  coefficients  at  these 
locations  were  0,  5,  3,  2,  1,  and  0,  respectively.  The  temperature 
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of  the  Plancklan  curve  Joining  the  last  two  knots  was  15E0. 

The  guessed  spectrum  had  six  knots  at  the  following  locations: 

0,  EO,  4E0,  35E0,  50E0,  and  128E0.  The  coefficients  at  these 
locations  were  0,  4,  4,  3,  2,  and  0,  respectively.  The  temperature 
of  the  Plancklan  curve  joining  the  last  two  knots  was  i3E0. 

2.  The  actual  and  guessed  cubic  splines  were  identical.  The 
spectra  have  the  same  knot  locations  and  coefficients  as  the 
actual  spectrum  of  test  case  1. 

The  cubic  spline  basis  function  method  produced  an  unfolded  spectrum, 
for  each  of  the  validation  cases,  that  was  nearly  identical  to  the  actual 
spectrum.  However  when  the  initial  guess  for  the  cubic  spline  basis 
function  was  chosen  a  little  less  carefully,  the  method  didn't  converge. 

Effects  of  Errors  on  the  Validation  Cases 

The  effect  of  measurement  and  calibration  errors  was  simulated  by  the 
addition  of  random  numbers  to  the  measured-to- predicted  ratio,  bi,  as  dis¬ 
cussed  in  the  previous  section.  A  demostration  of  this  effect  on  the  vali¬ 
dation  cases  is  illustrated  by  Figure  X.  In  this  figure,  the  modified 
iteration  method  was  run  in  the  presence  of  noise  for  case  2  (the  guessed 

spectrum  was  correct).  The  collection  of  unfolded  spectra  has  an  envelope 
which  follows  the  actual  curve,  but  each  unfolded  spectrum  has  details 

which  could  have  been  interpreted  as  being  significant. 
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The  dips  in  the  unfolded  spectra  happened  in  the  intervals  where  the 
noise  from  overlapping  detectors  were  both  negative.  The  peaks  were 
where  the  noise  from  the  overlapping  detectors  were  both  positive.  When 
the  noise  was  positive  for  one  detector  and  negative  for  the  overlapping 
detector,  the  overall  effect  of  the  noise  was  negligable. 

There  was  another  effect  of  the  noise  that  was  pecticular  to  the  Iter¬ 
ation  method.  The  iteration  method  was  run  using  a  geometric  spacing  of 
the  energy  bins  in  order  to  significantly  shorten  the  run  times.  For  the 

iteration  method  for  cases  with  no  noise,  the  difference  between  the  spec¬ 
tra  using  linear  spacing  and  geometric  spacing  was  negligable. 
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The  noisy  signal  produced,  for  the  geometric  energy  bins,  very  short 
spikes  for  small  energies  and  long  spikes  for  the  larger  energies.  This  is 
an  artifact  of  the  unfolding  method.  In  order  to  give  a  realistic  compari¬ 
son  for  the  noisy  signal  data  produced  by  the  iteration  method,  the  linear 
energy  bin  spacing  was  used. 

Test  Cases 

There  were  four  actual  spectra  used  to  compare  the  unfolding  methods. 
The  following  is  a  list  of  the  actual  spectra  used  for  the  test  cases: 

1.  Two  Planckians  with  a  spike.  This  spectrum  was  chosen  to 
compare  the  ability  of  the  methods  to  recover  broad  structure 
(the  two  Plankians)  with  their  ability  to  recover  finer  detail 
(the  spike). 

The  Planckians  were  T1  =  2E0,  coefficient  1  =  2,  T2  = 

10E0,  and  coefficient  2  =  5.  The  spike  was  11/8  EO  wide 
centered  at  10E0. 

The  guessed  spectrum  was  the  Planckian  spectrum  T  =  10E0 
and  coefficient  =  10.  For  the  Planckian  basis  function  method, 
a  second  two  temperature  Planckian  was  also  used. 

2.  Heaviside  function.  This  spectrum,  and  the  following  two, 
represent  the  ability  of  the  unfolding  method  to  represent 
discontinuous  behavior  of  the  actual  spectrum.  The  actual 
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spectrum  was  1  for  energies  less  than  10E0  and  0.5  for  energies 
greater  than  10E0. 

The  guessed  spectrum  for  all  methods  was  0.76  for  ail 
energies.  The  Plankian  basis  function  wasn’t  used  for  this  or 
the  following  two  discontinuous  actual  spectra. 

3.  Barrier  function.  The  actual  spectrum  was  1  between  energies 
8E0  and  12E0,  but  0.5E0  for  all  other  energies.  The  guessed 
spectrum,  for  all  the  methods,  was  1  for  all  energies. 

4.  Triangle  function.  The  actual  spectrum  was  an  isosceles 
triangle  with  an  apex  of  1  at  10E0  and  a  value  of  0.5  at  the 
base.  The  base  extended  from  6E0  to  14E0.  The  actual 
spectrum  was  0.5  for  all  other  energies.  The  guess  spectrum, 
for  all  the  methods,  was  1  for  all  energies. 

Both  the  modified  iteration  and  iteration  methods  made  an  attempt  to 
represent  the  spike  by  a  broader  shape  at  the  same  energy  as  shown  in 
Figures  XI  and  XII,  respectively.  The  pulse  seems  to  disrupt  the  spectrum 
for  larger  energies.  This  Is  probably  due  to  the  overlapping  response 
functions  in  the  energy  range  of  the  spike. 


FIGURE  XI:  Modified  Iteration,  Test  Case  1 
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The  Planckian  basis  function  method  produced  an  unfolded  spectrum 
composed  of  two  Plancklans  of  different  temperatures.  Because  of  the  rela¬ 
tive  temperatures  and  coefficients,  the  spectrum  appears  to  be  a  one  tem¬ 
perature  Planckian  as  shown  in  Figure  XIII. 
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FIGURE  XIII:  Planckian  Basis  Function,  Test  Case  1 
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The  two  Iteration  methods  have  a  ragged  approximation  to  the  Heavi 
side  spectrum,  rounding  off  the  corners  and  crossing  the  discontinuity  at 
the  midpoint,  Figure  XIV  and  XV.  For  the  cubic  spline,  the  unfolded 
spectrum  can’t  be  distinguished  from  the  actual  spectrum.  This  is  not 
unexpected  as  the  cubic  spline  has  a  good  ability  to  represent  straight 


lines. 
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For  test  case  3,  the  two  iteration  methods  move  their  unfolded  spectra 
to  follow  the  actual  spectrum  Into  the  spike  discontinuity.  Figure  XVI  and 
Figure  XVII.  The  cubic  spline's  unfolded  spectrum  with  seven  knots  follows 
the  discontinuous  spike  of  the  actual  spectrum  as  shown  in  Figure  XVIII. 
Except  for  the  dropoff  for  the  higher  energies,  the  cubic  spline  is  a  better 
approximation  than  the  iteration  methods. 


49 


FIGURE  XVI:  Modified  Iteration,  Test  Case  3 
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Both  iteration  methods,  follow  the  triangular  actual  spectrum  rounding 
off  the  apex  and  edges  of  the  trlan  -ie  similarly  to  the  representation  in 
case  3  as  seen  In  figures  XIX  and  XX.  For  the  discontinuities  studied,  the 
iteration  methods  are  not  that  particular  of  the  shape  of  the  spike.  The 
cubic  spline  Is  a  good  representation  of  the  actual  spectrum  (figure  XXI), 
but  it  still  has  the  high  energy  dropoff  seen  in  the  previous  case. 
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Noise  on  the  Test  Cases 

For  each  of  the  test  cases  and  unfolding  methods,  ten  unfolded  spectra 
were  generated  to  determine  the  bars  In  Figures  XXII  and  XXIII.  In  each 
of  the  figures  the  actual  spectrum  is  graphed  with  the  bars  superimposed 
upon  the  spectrum  to  represent  the  spread  of  spectra  due  to  noisy  signals. 

In  figure  XXII,  the  Planckian  method  has  smaller  spreads  in  the  noisy 
spectra  than  the  Iteration  method.  But  the  iteration  method  has  it's  most 
significant  spread  near  the  pulse  which  the  Planckian  method  almost 

Ignores. 
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In  figure  XXIII,  the  cubic  spline  and  iteration  methods  have  compara¬ 
ble  spread  in  their  noisy  spectra.  Both  methods  have  their  largest  spread 
near  the  spike,  which  both  attempt  to  represent,  and  decrease  with 
increasing  energy.  However  the  cubic  spline  still  has  large  variations  for 
large  energy,  where  the  spread  of  unfolded  spectra  from  the  iteration 
method  goes  to  nearly  zero. 
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V.  Summary  and  Conclusion 


Summary 

Validation  of  the  Methods.  The  purpose  of  the  validation  cases  is  to 
ensure  that  the  basis  function,  iteration,  and  the  modified  iteration  meth¬ 
ods  were  working  correctly.  Validation  cases  one  through  three  for  the 
basis  function  methods  produced  unfolded  spectra  which  were  nearly  an 
exact  match  to  the  actual  spectra. 

The  iteration  and  modified  iteration  methods  produced  spectra  which 
had  discontinuities  near  the  k-edges  of  the  response  function.  However 
the  overall  shape  was  a  good  representation  of  the  actual  spectra. 

Test  Cases.  In  case  one,  composed  of  two  temperature  Planckians  and 
a  spike,  the  Planckian  basis  method  proved  to  be  a  bad  representation  of 
the  actual  spectrum.  The  method  didn't  indicate  that  the  pulse  was  there 
and  the  overall  representation  was  that  of  nearly  a  one  temperature 
Planckian.  The  iteration  methods,  on  the  other  hand  proved  to  be  more 
adept  at  approximating  the  shape  of  the  spike.  For  case  four,  the  itera¬ 
tion  methods  and  the  cubic  spline  were  both  good  at  representing  the  dis¬ 
continuous  shape  of  the  spectrum.  The  unfolded  spectrum  of  the  cubic 
spline  was  better  at  representing  the  actual  form  the  the  spike,  but  it  had 
a  significant  dropoff  for  higher  energies. 

The  test  cases  for  the  noisy  measurement  signal  demonstrated  that  the 
basic  shape  of  the  unfolded  spectra  varies  considerably  from  trial  to  trial, 
even  if  there  is  no  change  in  the  data.  The  spread  in  the  spectra  is  very 
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small  for  the  Planckian  method,  but  it  doesn't  represent  the  discontinuity 
of  the  spike.  The  spread  for  the  cubic  spline  approximates  that  of  the 
iteration  methods  for  low  energies,  but  remains  large  for  large  energies. 

Conclusions: 

The  results  from  the  validation  and  test  cases  indicates  the  following 
conclusion  can  be  drawn: 

1.  Details  of  unfolded  spectra,  from  real  data  produced  during 
underground  nuclears  effects  tests,  can  not  be  reasonably 
assumed  to  reflect  real  details  in  the  actual  spectra.  Either 
the  details  of  the  actual  spectra  do  not  survive  the  unfolding 
process  in  the  presence  of  noise  or  the  details  on  the  unfolded 
spectra  are  added  during  the  unfolding  process. 

2.  The  actual  spectrum  can  not  be  unfolded  exactly,  even  if 
there  were  no  measurement  noise. 

3.  For  actual  spectra  which  is  very  nearly  Planckian,  the 
Planckian  basis  function  method  is  the  best.  However  for  the 
test  cases  in  this  study,  the  iteration  methods  do  better  than 
the  Planckian  method. 

4.  The  cubic  spline  method  is  comparable  to,  if  not  better  than, 

the  iteration  method  for  representing  discontinuous  spectra. 

It  does  have  significant  a  dropoff  for  higher  energies  and  a 

larger  spread  of  spectra  with  a  noisy  signal. 
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VI.  Recommendations 


The  following  recommendations  are  for  continuing  work  in  this  area. 
First,  determine  the  impact  of  more  detail  in  the  specification  of  the 
detectors  e.g.,  different  standard  deviations  and  response  functions. 

Secondly,  the  development  of  another  method  to  make  the  unfolded 
spectra  more  reasonable  looking.  Perhaps  a  scheme  where  the  local  points 
on  the  unfolded  spectra  are  balanced  to  make  the  behavior  of  the  large 
scale  less  abrupt. 

Thirdly,  the  use  of  basis  functions,  other  than  Plankian  and  cubic 
splines,  to  represent  the  spectra. 

Lastly,  determine  how  the  shape  and  size  of  the  pulse  perturbation 
can  affect  the  unfolded  spectrum  and  if  there  is  a  minimum  size  and  shape 
perturbation  which  can  be  unfolded. 
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Abstract 


The  purpose  of  this  study  is  to  ccrpare  methods  for  solving  the  Fred¬ 
holm  integral  equation  of  the  first  kind.  The  Fredholm  equation  has  several 
practical  applications  including  geology,  superconductivity ,  and  aerodynamics. 
Of  specific  interest  is  its  application  to  determining  radiation  spectra  using 
data  from  underground  nuclear  effects  simulations. 

The  two  basic  methods  studied  were  the  basis  function  method  and  the  iter¬ 
ation  method.  The  basis  function  method  is  a  representation  of  the  unfolded 
spectrum  by  a  series  of  Planckian  or  cubic  spline  functions.  The  iteration 
method  scales  the  unfolded  spectrum  so  that  its  weighted  integral  over  a  given 
interval  matches  that  of  the  actual  spectrum . 

Both  the  basis  function  methods  produced  excellent  results  when  the 
actual  spectrum  was  a  sum  of  it's  basis  functions.  The  cubic  spline  method 
produced  unfolded  spectra  which  were  good  approximations  for  discontinuous 
actual  spectra.  However,  there  was  a  significant  dropoff  of  the  spectrum 
for  the  cubic  spline  unfolded  spectra  for  higher  energies. 

The  iteration  methods  produced  accurate  approximations  for  actual  spectra 
that  were  both  basis  functions  and/or  discontinuous  spectra.  There  were  two 
problems  with  this  method:  the  iteration  method  produced  spectra  which  were 
discontinuous  at  the  discontinuities  of  the  weighting  function  and  noisy  data 
can  produce  large  discontinuities  in  the  unfolded  spectrum. 

For  the  cases  studied,  the  iteration  method  is  the  best  method,  because 
of  it's  ability  to  unfold  a  variety  of  spectra  (in  contrast  to  the  Planckian 
basis  function  method)  and  it’s  good  behavior  at  the  higher  energies  for  noisy 
data  (in  contrast  to  the  cubic  spline  method) . 
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